'%!in%' <- function(x,y)!('%in%'(x,y))

library(ggplot2)
library(stargazer)

# Load Data ---------------------------------------------------------------
load("lc_dat.rdata")

# Descriptives on Parenting Decisions -------------------------------------
# * Table 1 ---------------------------------------------------------------
# Talk  since May
sum(prop.table(table(lc_dat$tellkid_equal))[2:4]) # 0.9090065
sum(prop.table(table(lc_dat$tellkid_whiteworkharder))[2:4]) # 0.362117
sum(prop.table(table(lc_dat$tellkid_othersfacediscrim))[2:4]) #0.7892293
sum(prop.table(table(lc_dat$tellkid_blackrespectpolice))[2:4]) #0.4419684
sum(prop.table(table(lc_dat$tellkid_historicalfigs))[2:4]) #0.7994429
sum(prop.table(table(lc_dat$tellkid_whitediscrim))[2:4]) #0.5942433
sum(prop.table(table(lc_dat$tellkid_whiteprivilege))[2:4]) #0.4930362

stargazer(lc_dat[c("tellkid_equal", "tellkid_historicalfigs", "tellkid_othersfacediscrim", "tellkid_whitediscrim", "tellkid_whiteprivilege",  "tellkid_blackrespectpolice", "tellkid_whiteworkharder")], omit.summary.stat = c("p25", "p75", "sd", "max")) #sets up table and provides means and columns to input mode and pct responding yes, some amount

#get mode
table(lc_dat$tellkid_equal) #3
table(lc_dat$tellkid_historicalfigs) #2
table(lc_dat$tellkid_othersfacediscrim) #2
table(lc_dat$tellkid_whitediscrim) #0
table(lc_dat$tellkid_whiteprivilege) #0
table(lc_dat$tellkid_blackrespectpolice) #0
table(lc_dat$tellkid_whiteworkharder) #0

##how many parents said none of these?
lc_dat$tellkid_none[lc_dat$tellkid_equal == 0 & lc_dat$tellkid_historicalfigs == 0 & lc_dat$tellkid_othersfacediscrim == 0 & lc_dat$tellkid_whitediscrim == 0 & lc_dat$tellkid_whiteprivilege == 0 & lc_dat$tellkid_blackrespectpolice == 0 & lc_dat$tellkid_whiteworkharder == 0] <- 1
lc_dat$tellkid_none[lc_dat$tellkid_equal > 0 | lc_dat$tellkid_historicalfigs > 0 | lc_dat$tellkid_othersfacediscrim > 0 | lc_dat$tellkid_whitediscrim > 0 | lc_dat$tellkid_whiteprivilege > 0 | lc_dat$tellkid_blackrespectpolice > 0 | lc_dat$tellkid_whiteworkharder > 0] <- 0

table(lc_dat$tellkid_equal, lc_dat$tellkid_none)
table(lc_dat$tellkid_historicalfigs, lc_dat$tellkid_none)
table(lc_dat$tellkid_othersfacediscrim, lc_dat$tellkid_none)
table(lc_dat$tellkid_whitediscrim, lc_dat$tellkid_none)
table(lc_dat$tellkid_whiteprivilege, lc_dat$tellkid_none)
table(lc_dat$tellkid_blackrespectpolice, lc_dat$tellkid_none)
table(lc_dat$tellkid_whiteworkharder, lc_dat$tellkid_none)
prop.table(table(lc_dat$tellkid_none))


##Percentage who say others face discrim at least once and Black people should respect police for things to go better
table(lc_dat$tellkid_othersfacediscrim, lc_dat$tellkid_blackrespectpolice)
prop.table(table(lc_dat$tellkid_othersfacediscrim, lc_dat$tellkid_blackrespectpolice))

lc_dat$saidothersfacediscrim[lc_dat$tellkid_othersfacediscrim == 0] <- 0
lc_dat$saidothersfacediscrim[lc_dat$tellkid_othersfacediscrim > 0] <- 1
table(lc_dat$tellkid_othersfacediscrim, lc_dat$saidothersfacediscrim)

lc_dat$saidblackrespectpolice[lc_dat$tellkid_blackrespectpolice == 0] <- 0
lc_dat$saidblackrespectpolice[lc_dat$tellkid_blackrespectpolice > 0] <- 1
table(lc_dat$tellkid_blackrespectpolice, lc_dat$saidblackrespectpolice)

table(lc_dat$saidothersfacediscrim, lc_dat$saidblackrespectpolice)
prop.table(table(lc_dat$saidothersfacediscrim, lc_dat$saidblackrespectpolice))


##Percentage who talk about white privilege at least once and say whites work harder
table(lc_dat$tellkid_whiteprivilege, lc_dat$tellkid_whiteworkharder)
prop.table(table(lc_dat$tellkid_whiteprivilege, lc_dat$tellkid_whiteworkharder))

lc_dat$saidwhitepriv[lc_dat$tellkid_whiteprivilege == 0] <- 0
lc_dat$saidwhitepriv[lc_dat$tellkid_whiteprivilege > 0] <- 1
table(lc_dat$tellkid_whiteprivilege, lc_dat$saidwhitepriv)

lc_dat$saidwhitework[lc_dat$tellkid_whiteworkharder == 0] <- 0
lc_dat$saidwhitework[lc_dat$tellkid_whiteworkharder > 0] <- 1
table(lc_dat$tellkid_whiteworkharder, lc_dat$saidwhitework)

table(lc_dat$saidwhitepriv, lc_dat$saidwhitework)
prop.table(table(lc_dat$saidwhitepriv, lc_dat$saidwhitework))



# * Table 2 ---------------------------------------------------------------
##actions, nonpolitical: book battery, new art, child's school, anti-racism webinars
prop.table(table(lc_dat$movie_divchar)) #0.4076493 
prop.table(table(lc_dat$book_divhist)) #0.3193277  
prop.table(table(lc_dat$buy_divchar)) #0.3264925 
prop.table(table(lc_dat$book_divchar)) #  0.2782446  
prop.table(table(lc_dat$book_discrim)) # 0.2567694  
prop.table(table(lc_dat$blmsign)) # 0.233209 
prop.table(table(lc_dat$book_teachrace)) #0.1951447  
prop.table(table(lc_dat$communitymeeting)) #0.1968284 
prop.table(table(lc_dat$attendblm)) # 0.1902985 
prop.table(table(lc_dat$antiracism_wkid)) #0.1837687 
prop.table(table(lc_dat$changeschool)) #0.1772388 
prop.table(table(lc_dat$antiracism_parenting)) #0.136194 

stargazer(lc_dat[c("movie_divchar", "buy_divchar", "book_divhist",  "book_divchar", "book_discrim",  "blmsign", "book_teachrace", "communitymeeting", "attendblm", "antiracism_wkid", "changeschool", "antiracism_parenting")], omit.summary.stat = c("N", "p25", "p75", "sd",  "min", "max"))


##how many parents do none of these?
lc_dat$do_none[lc_dat$movie_divchar == 0 & lc_dat$book_divhist == 0 & lc_dat$buy_divchar == 0 & lc_dat$book_divchar == 0 & lc_dat$book_discrim == 0 & lc_dat$blmsign == 0 & lc_dat$book_teachrace == 0 & lc_dat$communitymeeting == 0 & lc_dat$attendblm == 0 & lc_dat$antiracism_wkid == 0 & lc_dat$changeschool == 0 & lc_dat$antiracism_parenting == 0 ] <- 1

lc_dat$do_none[lc_dat$movie_divchar == 1 | lc_dat$book_divhist == 1 | lc_dat$buy_divchar == 1 | lc_dat$book_divchar == 1 | lc_dat$book_discrim == 1 | lc_dat$blmsign == 1 | lc_dat$book_teachrace == 1 | lc_dat$communitymeeting == 1 | lc_dat$attendblm == 1 | lc_dat$antiracism_wkid == 1 | lc_dat$changeschool == 1 | lc_dat$antiracism_parenting == 1 ] <- 0

table(lc_dat$movie_divchar, lc_dat$do_none)
table(lc_dat$buy_divchar, lc_dat$do_none) 
table(lc_dat$book_divhist, lc_dat$do_none)
table(lc_dat$book_divchar, lc_dat$do_none)
table(lc_dat$book_discrim, lc_dat$do_none)
table(lc_dat$blmsign, lc_dat$do_none)
table(lc_dat$book_teachrace, lc_dat$do_none)
table(lc_dat$communitymeeting, lc_dat$do_none)
table(lc_dat$attendblm, lc_dat$do_none)
table(lc_dat$antiracism_wkid, lc_dat$do_none)
table(lc_dat$changeschool, lc_dat$do_none) 
table(lc_dat$antiracism_parenting, lc_dat$do_none)
prop.table(table(lc_dat$do_none))

##number of things done: 
lc_dat$numacts <- lc_dat$movie_divchar + lc_dat$buy_divchar + lc_dat$book_divhist + lc_dat$book_divchar + lc_dat$book_discrim + lc_dat$blmsign + lc_dat$book_teachrace + lc_dat$communitymeeting + lc_dat$attendblm + lc_dat$antiracism_wkid + lc_dat$changeschool + lc_dat$antiracism_parenting
table(lc_dat$numacts)
prop.table(table(lc_dat$numacts))
sum(prop.table(table(lc_dat$numacts))[1]) #0 acts: 29.63%
sum(prop.table(table(lc_dat$numacts))[2:13]) #at least 1 act: 70.37%
sum(prop.table(table(lc_dat$numacts))[3:13]) #at least 2 acts: 56.45%
sum(prop.table(table(lc_dat$numacts))[4:13]) #at least 3 acts: 46.17%
sum(prop.table(table(lc_dat$numacts))[5:13]) #at least 4 acts: 32.15%
sum(prop.table(table(lc_dat$numacts))[6:13]) #at least 5 acts: 27.01%
sum(prop.table(table(lc_dat$numacts))[7:13]) #at least 6 acts: 21.40%
sum(prop.table(table(lc_dat$numacts))[8:13]) #at least 7 acts: 15.05%
sum(prop.table(table(lc_dat$numacts))[9:13]) #at least 8 acts: 10.47%


# Correlates of Parenting Decisions Models --------------------------------
# * Specify DVs ------------------------------------------------------------
# combining to DV list
DVS <- c("book_divchar", "book_divhist", "book_discrim",              
         "book_teachrace", "tellkid_equal", "tellkid_whiteworkharder",
         "tellkid_othersfacediscrim", "tellkid_blackrespectpolice", "tellkid_historicalfigs",
         "tellkid_whitediscrim", "tellkid_whiteprivilege", "buy_divchar",
         "movie_divchar", "changeschool", "attendblm",                 
         "blmsign", "communitymeeting", "antiracism_wkid",
         "antiracism_parenting")

# * Specify RHS -------------------------------------------------------------
RHS <- c("age_lucid_sc", "woman", "edu_sc", "inc_sc",
         "zip_white_prop", "zip_black_prop", "zip_med_inc_sc", "zip_college_prop",
         "oldest_child_age_sc", "multkids", "childcaretime_sc", "sahp", 
         "covidwork_hours", "covidwork_left",
         "pid7_sc", "ideo7_sc")


# ** Model Loop --------------------------------------------------------------
# Formatting DF to store coefficients and p-values
coef_mat <- array(NA, c(length(DVS), length(RHS) + 1))
coef_mat <- as.data.frame(coef_mat)
names(coef_mat) <- c("var", RHS)
coef_mat$var <- DVS
sig_mat <- coef_mat

# empty list to store individual model results
model_results <- list()

for(i in 1:nrow(coef_mat)){
  rhs <- paste(RHS, collapse = "+")
  m <- lm(as.formula(paste(coef_mat$var[i], rhs, sep = "~")), 
          data = lc_dat)
  m_coef <- summary(m)$coefficients
  # store coefficients excluding intercept
  coef_mat[i, -1] <- m_coef[-1, 1]
  # store p-values but for intercept
  sig_mat[i, -1] <- m_coef[-1, 4]
  
  model_results[[i]] <- m
}

# *** Format DF to Plot -----------------------------------------------------
mod_out <- coef_mat
mod_out <- reshape::melt(mod_out, id = c("var"))
names(mod_out)[3] <- "coef"
# reshaping to merge
sig_mat2 <- reshape::melt(sig_mat, id = c("var"))
names(sig_mat2)[3] <- "p_value"
mod_out <- merge(mod_out, sig_mat2, by = c("var", "variable"))

# for plotting
mod_out$sign <- ifelse(mod_out$coef < 0, "Negative", "Positive")
mod_out$sig <- ifelse(mod_out$p_value <= 0.05, "*", "")

# ** Figure 3 ----------------------------------------------------------------
iv_labs <- c("Age",
             "Woman",
             "Education",
             "Income",
             "Zip: Prop White",
             "Zip: Prop Black",
             "Zip: Median Income",
             "Zip: Prop College Degree",
             "Oldest Child's Age",
             "Multiple Kids",
             "Childcare Time",
             "Stay at home Parent",
             "COVID Reduced Hours",
             "COVID Left Workforce",
             "Partisanship (Republican)",
             "Ideology (Conservative)")

dv_labs <- c("Anti-Racist\nParenting Workshop",
             "Anti-Racist\nWorkshop w/kid",
             "Attend BLM Protest",
             "Make BLM Sign",
             "Book on\nDiscrimination",
             "Book w/\ndiverse characters",
             "Historical Books",
             "Book on\nTeaching about Race",
             "Toys w/diverse characters",
             "Change School",
             "Attend\nCommunity Meeting",
             "Movie w/diverse characters",
             "Tell Kids: Respect Police",
             "Tell Kids: People Equal",
             "Tell Kids:\nHistorical Figures",
             "Tell Kids:\nOthers Face Discrimination",
             "Tell Kids:\nWe Face Discrimination",
             "Tell Kids:\nWhite Privilege",
             "Tell Kids:\nWhites work Harder")

mod_out$dv <- factor(mod_out$var,
                     labels = dv_labs)
mod_out$dv <- factor(mod_out$dv,
                     levels = levels(mod_out$dv)[c(14:16, 13, 19, 17, 18, 6, 9, 12, 5, 7, 1, 2, 8, 10, 3, 4, 11)])
# table(mod_out$dv, mod_out$var)


mod_out$iv <- mod_out$variable
levels(mod_out$iv) <- iv_labs
# table(mod_out$iv, mod_out$variable)

mod_out$abs_coef <- abs(mod_out$coef)
mod_out$sign_symb <- ifelse(mod_out$sign == "Negative", "-", "+")
mod_out$pos <- ifelse(mod_out$sign == "Positive", "Positive", NA)
hw <- .9

# figure
ggplot(mod_out, aes(x = dv, y = iv, fill = coef)) +
  geom_tile(aes(color = neg), width = hw, height = hw, size = 1) +
  scale_fill_gradient2(low = "darkorange1", 
                       high = "navy",
                       name = "Coefficient\nEstimate",
                       guide = "legend",
                       n.breaks = 6, limits = c(-.55, .55)) +
  scale_color_manual(name = "", values = "black", guide = FALSE) +
  geom_text(aes(x = dv, y = iv, label = sig), 
            color = "black", size = 5) +
  ggridges::theme_ridges() +
  labs(x = "", y = "",
       caption = "* p < .05. Boxes highlight negative coefficient estimates.") +
  scale_x_discrete(position = "top") +
  scale_y_discrete(limits = rev) +
  theme(axis.text.x = element_text(angle = 45, hjust = .1, size = 11,
                                   lineheight = .65),
        legend.position = "bottom",
        legend.text = element_text(size = 12),
        plot.margin=unit(c(0,2,0,0),"cm"))
# ggsave("fig3.pdf", width = 10, height = 9)

# ** Appendix D Tables -------------------------------------------------------------------
iv_labs <- c("Age",
             "Woman",
             "Education",
             "Income",
             "Zip: Prop White",
             "Zip: Prop Black",
             "Zip: Median Income",
             "Zip: Prop College Degree",
             "Oldest Child's Age",
             "Multiple Kids",
             "Childcare Time",
             "Stay at home Parent",
             "COVID Reduced Hours",
             "COVID Left Workforce",
             "Partisanship (Republican)",
             "Ideology (Conservative)")
dv_labs <- c("Anti-Racist Parenting Workshop",
             "Anti-Racist Workshop w/kid",
             "Attend BLM Protest",
             "Make BLM Sign",
             "Book on Discrimination",
             "Book w/ diverse characters",
             "Historical Books",
             "Book on Teaching about Race",
             "Toys w/diverse characters",
             "Change School",
             "Attend Community Meeting",
             "Movie w/diverse characters",
             "Tell Kids: Respect Police",
             "Tell Kids: People Equal",
             "Tell Kids: Historical Figures",
             "Tell Kids: Others Face Discrimination",
             "Tell Kids: We Face Discrimination",
             "Tell Kids: White Privilege",
             "Tell Kids: Whites work Harder")

# model results by DV order
stargazer(model_results[c(19,18,15,16,3,1)],
          title = "Correlates of Parenting Decisions",
          dep.var.caption = "",
          covariate.labels = iv_labs,
          column.labels = dv_labs[1:6],
          dep.var.labels.include = F,
          header = F, initial.zero = F,
          type = "latex", label = "heat_map_tab1=", 
          out = "table1.tex",
          digits = 3, omit.stat = c("f", "adj.rsq"),  df = F, no.space = T)

stargazer(model_results[c(2, 4, 12, 14, 17, 13, 8)],
          title = "Correlates of Parenting Decisions",
          dep.var.caption = "",
          covariate.labels = iv_labs,
          column.labels = dv_labs[7:13],
          dep.var.labels.include = F,
          header = F, initial.zero = F,
          type = "latex", label = "heat_map_tab2", 
          out = "table2.tex",
          digits = 3, omit.stat = c("f", "adj.rsq"),  df = F, no.space = T)

stargazer(model_results[c(5, 9, 7, 10, 11, 6)],
          title = "Correlates of Parenting Decisions",
          dep.var.caption = "",
          covariate.labels = iv_labs,
          column.labels = dv_labs[14:19],
          dep.var.labels.include = F,
          header = F, initial.zero = F,
          type = "latex", label = "heat_map_tab3", 
          out = "heat_map_table3.tex",
          digits = 3, omit.stat = c("f", "adj.rsq"),  df = F, no.space = T)
